Beam oscillations and curling in chirped periodic structures with metamaterials 
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We study the propagation of electromagnetic waves in one-dimensional chirped periodic structures 
composed of alternating layers of negative-index (or left-handed) metamaterial and conventional 
dielectric, under the condition of the zero average refractive index. We consider the case when the 
periodic structure has a linear chirp, and the chirp is introduced by varying linearly the thickness 
of the layers across the structure. We apply an asymptotic analytical method for the analysis of the 
Bloch oscillations and find that, in a sharp contrast to ordinary periodic dielectric structures, the 
energy flow in multi-layered stacks with metamaterials may have the opposite direction at the band 
edges, thus providing novel possibilities for the beam steering in the transmission band. 



I. INTRODUCTION 

Materials with negative index of refraction, also known 
as left-handed metamaterials, are attracting a great sci- 
entific interest nowadays. This novel type of materials 
was predicted theoretically back to 1967 [l[ but only 30 
years later this theoretical curiosity was followed by the 
first experimental demonstrations based on the fabrica- 
tion of composite structures consisting of split-ring res- 
onators and metallic wires 0, H|- The basic properties 
of the left-handed metamaterials predicted theoretically 
have been confirmed by experiment, and this led to a 
rapid progress in this field and many other interesting 
discoveries. In particular, many new physical phenomena 
based on the concept of metamaterials and negative in- 
dex of refraction have been predicted, including the fun- 
damental concept of perfect lens [i| and, more recently, 
electromagnetic cloaking @, H, 0] • 

One of the interesting directions in the analysis of the 
unusual properties of this novel composite metamateri- 
als is the study of periodic structures (or photonic crys- 
tals) composed of metamaterials. Photonic crystals, be- 
ing widely used for both light control and beam manip- 
ulation, demonstrate a variety of novel physical effects, 
and they may offer many new possibilities being com- 
bined with metamaterials. In particular, a series of recent 
studies revealed the existence of a novel zero-index ban- 
gaps [1, [l(| [IH associated with one-dimensional pho- 
tonic crystals composed of alternating layers of negative- 
index and conventional dielectrics for which the aver- 
aged refractive index vanishes. Such novel periodic struc- 
tures demonstrate many intriguing properties, including 
substantial suppression of the Anderson localization and 
long- wavelength resonances ■ 

In this paper, we study the propagation of electromag- 
netic weaves in one-dimensional periodic structures com- 
posed of alternating layers of left-handed metamaterial 
and conventional dielectric. It is well established [HI, [l4| 
that introducing a chirp into a one-dimensional periodic 
structure leads to the generation of optical Bloch oscil- 
lations which can be also predicted in the presence of 
metamaterials [l5j . In addition to our numerical analysis 
of the Bloch oscillations studied earlier [Tj|, in this pa- 



per we demonstrate a number of novel, unique features 
of the beam propagation in this kind of chirped peri- 
odic composite structures. In particular, in contrast to 
Ref. [l5|, here we consider the chirped structures with a 
large number of layers and slowly varying period. Apply- 
ing the approximation based on the geometric optics, we 
reveal a variety of novel effects that could be observed in 
such structures, including the reversal of the energy flow 
at the band edges and the fascinating effect of the beam 
curling. This novel effect of the beam curling can be use- 
ful for the beam steering in the transmission band. We 
also study the eigenvalue problem and the corresponding 
eigenmodes and determine both the field structure and 
Poynting vector distribution associated with the beam 
propagation. 

The paper is organized as follows. In Sec. [IT] we intro- 
duce our structure as a one-dimensional stack of layers 
composed of alternating layers of metamaterial and di- 
electrics. We discuss the way of presenting eigenmodes 
of the structure with a slowly varying period and discuss 
the initial conditions. In Sec. lIIII we implement the Bloch 
theorem and derive the dispersion relation for the Bloch 
wave vector. Then, we analyze the bandgap diagram for 
this structure and discuss the conditions for the obser- 
vation of a novel type of the beam dynamics associated 
with the Bloch oscillations. Section HVl is devoted to the 
geometric optics approximation. Here we demonstrate 
a possibility to apply the geometric optics for the anal- 
ysis of layered structures with slowly changing parame- 
ters. We derive the equations of motion for the paraxial 
beams, and also study the period of oscillations pointing 
out novel properties of the beam trajectories connected 
with a change of the sign of the spatial group velocity 
in the structure due to the presence of a metamaterial. 
Also, we analyze different regimes of the beam propaga- 
tion and compare them with the well-known case corre- 
sponding to the similar structures composed of conven- 
tional dielectric layers. In Sec. [V] we discuss the eigen- 
value problem and point out the difficulties and failure 
of the commonly used numerical techniques. In addition, 
we describe an alternative method to find the eigenmodes 
of the structures with a large number of periods. Based 
on this method, we study the field distribution in the 
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structure and determine the direction of the Poynting 
vector in the beam. We observe a close correspondence 
with predictions of the geometric optics. Finally, Sec. VI 
concludes the paper. 



II. STRUCTURE GEOMETRY AND 
EIGENMODE FORMALISM 

We study a one-dimensional chirped periodic struc- 
ture shown schematically in Fig. [TJ where the slabs of 
a negative-index metamaterial with the width dim are 
separated by the layers of the conventional dielectric with 
the width d r m. We describe the variation of the refrac- 
tive index in the m-th pair of layers as follows, 



n{z) = 



m = - Jem 



where ni and n r are the refractive indices of metamaterial 
and dielectric layers, respectively, A m is the width of the 
m-th unit cell. In contrast to our previous studies [Til ], 
here we analyze the structures with a large number of 
unit cells and slowly varying period. 

We consider TE-polarized waves with the electric field 
described by one component E = (£^,0,0), and the 
waves propagating in the (y, z) plane. In this case the 
problem is described by the Helmholtz equation, 



A 2 E x {y,z)+n(zYE x (y,z) 



1 dn(z) dE x (y,z) 



dz 



dz 



= 0, (2) 



where A2 is the two-dimensional Laplacian. The field is 
assumed to be monochromatic with the frequency uj, and 
the coordinates are expressed in the units of c/uj, c is the 
speed of light. 

For the structures with slowly varying periods, the 
cigcnmodcs of the problem @ can be introduced by ap- 
plying the Bloch theorem for periodic systems and the 
corresponding Bloch- wave formalism [161 ] . In this ap- 
proach, the electric field corresponding to the eigenmode 
with the propagation constant k y in the unit cell with 
the width A can be represented as, 



E A (z, y) = A A U A (z) exp [-i{K^z + k y y)] , 



(3) 



where k y and Kb are normalized to lu/c, the amplitude 
Aa varies slow across the structure, U\ is the Bloch func- 
tion, and is the Bloch number of the periodic struc- 
ture with the unit cell size A, which we analyze system- 
atically in Sec. IIIII 

We study propagation of the Gaussian beams with the 
electric field at y = defined as 

E x (z, y = 0) = cxp[-(z - z j 2 /h]E(z) = %l>(z)E{z), (4) 



where E(z) is the eigenmode of the Eq. ([3]) with the prop- 
agation constant k y = ky, which will determine propaga- 
tion direction; h is the beam width, which we consider 
to be significantly large, i.e. h 3> < A >. The field in 
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FIG. 1: Top: Schematic of the chirped layered structure with 
a linearly growing period A m . Shaded slabs correspond to 
metamaterial layers with the width d; m separated by dielec- 
tric layers with the width d rm - Bottom: Schematic of the 
period variation across the structure. Structure with a lin- 
early changing period is embedded into semi-infinite periodic 
structures with the periods A; and A r , respectively. 



the structure can be represented as a superposition of 
eigenmodes. As long as the Gaussian beam is wide with 
respect to the unit cell, the spectrum of the contributing 
eigenmodes is narrow, and it is centered near the point 
ky [l6| • Such beam can be treated within the paraxial ap- 
proximation. The beam propagation is mainly described 
by the modes with k y ~ k y and the tangent of the beam 
angle is defined by (dk y /dK b )\ ky=k o A=Ao , where A cor- 
responds to the beam launching point. An additional 
beam tilt can be introduced in more general case of com- 
plex ip{ z )- 



III. PERIODIC STRUCTURE: BANDGAP 
PROPERTIES 



As the first step of our analysis, we study the periodic 
structure with the period A without a chirp. In this case 
Eq. ([2]) has periodic coefficients, and we can apply the 
Bloch theorem [l7| and present the electric field in the 
form, 



E x (z,y) = U(z) exp[-i(K b z + kyy)], 



(5) 



where Kb is the Bloch number, and the field envelope 
U(z) is a periodic function with the period A, so that 
U(z + A) = U(z). In the m-th pair of layers, the Bloch 
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FIG. 2: (Color online) Left: Bloch-wave phase incursion in 
the transmission regime for different values of the propagation 
constant k y . Inset shows a magnified part of the plot with the 
mode crossing. Middle: bandgap diagram for e r = fj. r = 1, 
a/b — 2 and ei = —5, fii = —0.8; low-order transmission reso- 
nance is shown. Right: group velocity dk y /dK^ for different 
values of the propagation constant. 

waves are presented as, 

UiAz) = [a^ r exp {-ik zLzr (z - mA)) + (6) 
bi <r exp (ik z i tZr (z — mA))] exp(zitb(;z — mA)), 

where the indices I and r correspond to the left- and right- 
handed slabs, respectively, and the amplitudes a^ r and 
bi tr are found from the boundary conditions at the inter- 
faces separatin g th e metamaterial and dielectric layers 
(see, e.g., Refs. [16l fl7l liljj). The Bloch wavenumber Kb 
is defined from the dispersion relation, 

2cos(AT;,A) = 2cos(k zr d r )cos(k z idi) — (7) 
fkzijM_ + kzrHi \ s ^ kzrdr ^ sin^d,), (8) 

\KzrfJLt Kzlfh-J 

where k z i }Zr = T J nf ~ — ky and k y is the propagation 
constant along the y axis. 

According to Eq. ([5|), the waves can propagate in the 
structure when Kb is real. In Fig. ^middle), we plot the 
bandgap diagram of the layered structure on the plane 
(A,k y ). Here we assume that the dielectric layer is air, 
e r = /i r = 1, and the dielectric layers are two times 
thicker than the metamaterial layers, i.e. d rm /di m = 
2, and this ratio is preserved across the structure. We 
choose the metamaterial with the parameters ei = — 5 
and /Lt; = —0.8. 

The bandgap properties of one-dimensional periodic 
structures with metamaterials have been studied in sev- 
eral papers [HI , including the case of the zero-index 
averaged refractive index |9j. It was shown that for 
the structures with J n(z)dz = the bandgap spectrum 

A 

includes transmission resonances which shrink into in- 
finitesimally thin lines into a complete bandgap when 



the widths of the slabs coincide [9(. For the normal in- 
cidence (i.e. when k y = 0), the transmission is observed 
only when n r d r = nidi = irq, where q is integer. 

We study the Bloch wave phase incursion Kf,A across 
the structure for different values of the propagation con- 
stant k y , see Fig. [SJleft). For small values of the propa- 
gation constant (e.g., k y < 0.5), the low-order transmis- 
sion resonances have zero phase incursion, in contrast to 
the case of conventional dielectric Bragg gratings. The 
maximum of the phase is accumulated in the middle of 
the transmission band, and its amplitude increases with 
the increase of the band order, reaching the value of 
7r in the n-th order transmission resonance. A growth 
of the propagation constant leads to the band coupling 
(0.5 < k y < 1), not shown here. The phase incursion in 
a new coupled broadband region is approaching ir. Fur- 
ther increase of the propagation constant leads to the 
narrowing of the transmission regions and to a rapid in- 
crease of the phase across the band to the value of ir. 
On the bandgap diagram, we observe two turning points 
at k y ~ 1.7 and k y ~ 2, (marked points 1 and 2 in the 
inset of Fig. ^middle)). At these points the modes with 
close k y have the same value of the Bloch wavenumber 
near the band edge. We trace the behavior of the band 
edges C_(fc y ) and C+(^y) with a change of the propa- 
gation constant k y . For k y ~ 1.7, the left band edge 
£_ changes the type of its monotonicity. Consequently, 
in this region the left (£_) and right (£+) band edges 
have different slopes, i.e. (c?C- / dk y ){dQ+ / dk y ) < 0, as 
shown in the inset of Fig. [3l The latter, with considera- 
tion that phase is growing from to n in the transmis- 
sion region, leads to the crossing of the corresponding 
modes (see inset in Fig. Finally for k y ~ 2, the 
right band edge C+ also changes the sign of its slope, and 
(d(- /dky)(d(+/dk y ) becomes positive. We note that for 
small values of k y , the first-order band edges also have 
different signs of their slopes but a new type of behavior 
cannot be observed because the phase incursion vanishes 
in the band for k y < 0.5. Also, in the conventional Bragg 
gratings both band boundaries have the same slope and 
thus the mode crossing does not occur. 

Crossing of two modes means that the effective spatial 
velocity defined as v = dkyjdK^., has different signs at 
different sides of the crossing point, and it becomes sin- 
gular at the point itself. As long as the velocity is related 
to the energy flow, we observe different directions of the 
energy flow in the y-direction at different edges of the 
band, as will be discussed in more details in Sec. II VI and 
Sec. |V] We plot the velocity for different values of the 
propagation constant k y , see Figs. [5Ja-c). For k y out- 
side the region, Fig. [2a), the velocity does not change 
its sign in the band, thus meaning that energy flow has 
the same direction at any point of the band. Similar 
type of the beam evolution can be found in the stacks 
of conventional dielectric layers. Further increase of k y 
leads to more asymmetric velocity profile and, finally, at 
k y ~ 1.7 [see Fig. EJb)] , the velocity becomes infinite at 
the left band edge. In Sec. IIVI below we demonstrate 
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that this case corresponds to the energy flow vanishing 
in the y-direction. When the propagation constant k y 
is inside the region [see Fig. Etc)], the velocity profile 
has the second-order discontinuity with the sign changed 
inside the band. Consequently, the energy flow in the 
y-direction also changes its sign within the transmission 
band. The singularity point in the spatial velocity pro- 
file moves from the left band edge to the right edge with 
the growth of the propagation constant from k y ~ 1.7 to 
k y ~ 2, (see turning points 1 and 2 in Fig. [3]). Finally, 
for the propagation constants k y > 2, the velocity pro- 
file becomes continuous and negative, indicating that the 
waves are backward in the whole band. 

Below, wc consider the simplest case of a linear chirp 
of the structure, Aj+i = A.; + SA with SA <C A, see Fig.[T] 
For 5A — ► 0, i.e. an adiabatic change of the unit cell size 
across the structure, we implement the geometric optics 
approximation. 



As long as the width of the unit cell changes adiabat- 
ically across the structure, the average velocity of the 
energy flow v e changes quasi-continuously. Therefore, 
we can define the rays in the space (A, y) (we note here 
that due to the linear chirp of the structure, the spatial 
coordinate z m is a linear function of A) as trajectories 
along which the average energy < S > is directed. Thus, 
the ray equation is defined as follows: 



dr 

dt 



(12) 



where r = e y y + e z A is the radius-vector describing the 
ray. 

In Ref. [l6[ it was shown that for a paraxial beam the 
average velocity of the energy flow coincides with the 
group velocity defined as 



duj 

= V g =VKO,= M , 



(13) 



IV. GEOMETRIC OPTICS APPROXIMATION 

We study the evolution of the Poynting vector by em- 
ploying an analogy with a homogeneous medium having 
gradually changing refraction index [l8|, [l!|. We con- 
sider the time-averaged Poynting vector, 



S = — Rc[E x H* 

07T 



(9) 



87T/1 



-Re i 



&zEx 



dE x 



d- 



dE x 
dy 



where E x is given by Eqs. ((U, e y and e z are the unit vec- 
tors. The time-averaged electromagnetic energy density 
is defined as 



-Lrc^ee* 

l07T 



Re sfi\E x 



dE* 



M HH*) 



dE x 



(10) 



dz 



dy 



In a homogeneous medium, the time-averaged Poynt- 
ing vector is tangential to the beam trajectory. It changes 
gradually with continuously varying refractive index, 
thus defining a geometric ray. In our case, the Poynt- 
ing vector changes its direction stepwise within each unit 
cell, since the refractive indices of the slabs comprising 
the unit cell have the opposite signs. To describe an 
average behavior of the energy flow on larger scale, we 
consider the Poynting vector and energy density averaged 

A 

over the unit cell, < S >= J Sdz and energy density 

o 

A 

< w >= 4- J wdz, respectively, 
o 

We introduce the average velocity of the energy 
flow USES: 



< S > 

< w > 



(11) 



where K = e y k y + e z Kb is the wave vector, and uj is the 
angular frequency. 

For a fixed frequency u> of monochromatic beam we ob- 
tain the dispersion relation k y = f(Kb), (see Eq. [7|); cor- 
responding to the curve in the space (k y ,Kb) described 
by the wave vector K, which is analogous to the equifrc- 
quency surface for anisotropic media [201 ] - Consequently, 
using Eqs. (fTTj) , (TT2"]) . and (flB")) we obtain v g = v e and 
< S > are normal to the curve k y = f(Kb). The latter 
condition is described by the equation, 



< S y > 
<S Z > 



dk y 



(14) 



where < S y > and < S z > are the y and z components 
of the Poynting vector < S > , v is the spatial velocity of 
the ray. 

Relation (fl"4"|) between the Poynting vector components 
and velocity v shows that when v — > oo, < S y >=0, and 
the total energy flows along the z-dircction. On the other 
hand, when v = 0, the energy flows along the structure 
(< S z >= 0). 

Using Eqs. ([Til) and OHi we rewrite Eq. (12]) as fol- 
lows: 



_ < S y > 
dt < w > 
dA _ <S Z > 
~dt ~ 



(15) 



< Sy > V 



<W> 



< w > 



Now we can derive the ray trajectory equation in the 
space (A, y) 



dA 
dy 



dk v 



dK h ' 



(16) 



where =F correspond to the forward and backward prop- 
agating waves, respectively. 

We note that the motion equation (fTo| can also be 
derived from the Format principle. We multiply Eq. |JJ 
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FIG. 3: (Color online) Left: period of beam oscillations. In- 
set shows the region of the sign change. Middle: bandgap 
diagram. Inset shows the turning points. Right: spectrum of 
the beam centered at the point ky = 1.843 



by K and using 
in the form [2011 



Eq. (j 1 3(1 . wc write the Format principle 



6 / Kdr = 0, 



(17) 



where dr — e y dy + e z dA is the unit vector along the 
beam trajectory. The Fermat principle is the principle 
of least action with the coordinate y playing the role of 
time. The Lagrange function L = ±Kb(dA/dy) + k y is 
defined accurate to the full y-derivative of A. As a result, 
the unambiguity in the Bloch wavenumbcr. Kb = K^ + 
27rm/A does not affect the ray behavior. The Hamilton 
equations for the rays can be derived in the form 19:]: 



dA 
dy 



dk y 



dK b 
dy 



±- 



dk y 
~dA 



(18) 



We remind now t licit ky remains constant across the 
structure, consequently dk y /dA = 0, meaning that the 
problem is invariant to translation along the layers. 
Hence, we consider only the first equation, which can be 
also derived from general suggestions, assuming that ky 
refers to the " energy" , being preserved in the structure 
and dky /dKb being the "group velocity" of the paraxial 
beam. 

The first equation of the system (jT5J) describes diffrac- 
tionlcss motion of a beam with a narrow spectrum cen- 
tered near the point ky. The beam motion is equivalent 
to the motion of an effective particle in a one-dimensional 
potential W = [v 2 {A)/2] between left (£_) and right ((+) 
band edges. As long as the spectrum of allowed energies 
near the point ky is quasi-continuous for infinite struc- 
tures with adiabatically changing width A, the motion 
remains periodic with the period of oscillations <f> found 
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FIG. 4: (Color online) Ray trajectories calculated in the ge- 
ometric optics approximation. Left: conventional dielectric 
structures, for /j, a = £a — 1 and = 1; e = 10. (a-e) Struc- 
tures with metamaterials for ky = 1.7, 1.85, 1.895, 1.95, 2, 
respectively. 



as. 



*(*v) = 2 



dA 



C- 



v{A, k y ) ' 



(19) 



where £_ and (+ refer to the band edges, v is defined by 
Eq. (fl"4|) . In the gap regions, the field decays exponen- 
tially across the structure so that the motion is prohib- 
ited, and the spatial velocity vanishes at the band edges. 
Consequently, the expression under the integral has sin- 
gular points at £±, and the integration fails. To resolve 
this problem, we introduce a new variable £, 



A 



C-+C+ , C--C+ 



and obtain 



<h 



cos(£), 



(C+-C-) 

2 KW 



(20) 



(21) 



The expression under the integral is finite within the inte- 
gration boundaries, hence the new integral is converging. 
Using the integral (|2"Tj) we calculate the dependence of the 
period on the propagation constant k y , see Fig. [3jlcft). 

We observe that with a growth of k y the oscillation 
period grows as well, this is explained by broadening of 
the band described in Sec. |TTJ A further increase of k y 
leads to the beam narrowing and decrease of the oscilla- 
tion period For 1.7 < k y < 2, i.e between the turning 
points, the period vanishes meaning that the correspond- 
ing mode does not transfer energy along the structure. A 
further increase of k y leads to negative values for the pe- 
riod of oscillations, which corresponds to the total energy 
flow in the direction opposite to the propagation constant 

Using Eq. (fl8|) and Eq. (|2"Tj) , we calculate the beam 
trajectories. First of all, we verify our results for the 
well-known structures with conventional dielectric slabs, 
and the trajectory corresponding to this case is shown in 
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Fig. rjjlcft). Wc assume that a dielectric layer with the 
width b has the parameters = 10 and jug, = 1, and 
the layers are separated by vacuum of the width a, so 
that a/b = 2. We reveal a close correspondence between 
our results and the trajectories calculated by means of 
the geometric optics in Ref. The beam is reflected 

from both band boundaries, where the spatial velocity 
vanishes. The total energy flow along the y-direction 
is positive, and it preserves its sign during oscillations. 
For the structures with metamaterials, the beam with 
the spectrum centered near the propagation constants 
k y outside of the region between turning points (points 
1 and 2 in Fig. [3]) experiences the same reflection from 
the boundaries of the Brillouin zone with the total energy 
flow along the structure in the positive direction, for k y < 
1.7, and in opposite (negative) direction, for k y > 2. 



Next, we study the beam propagation for the values of 
k y inside the region between the turning points. The tra- 
jectory calculated near the first turning point (see Fig. [4]) 
has a peculiarity near the left band edge. Correspond- 
ing period of oscillations is shown in Fig. [3^1eft) . Near 
the left band edge the group velocity becomes infinite, 
thus energy flow along the layers vanishes, and the energy 
flow across the structure changes its sign. Hence a vor- 
tex structure is formed at this point. The observed spike 
is formed by narrowing of the trajectory near the left 
band edge with increase of propagation constant until it 
reaches the first turning point. Further increase of prop- 
agation constant, as was described in Sec. [ill corresponds 
to a shift of the singular point of the group velocity fur- 
ther into the band region. The trajectories corresponding 
to this case are shown in Figs.BJa-e). Since the group ve- 
locity changes its sign before and after the singular point, 
the trajectory crosses itself and a beam curl is formed. 
For ky = 1.85 wc observe a curl near the left band edge 
of the structure, see Fig. 0Jb). Further increase of the 
propagation constant leads to the increase of the curl's 
size, and decrease of the period, see Fig. [3] Finally, at 
ky ~ 1.895 the curl reaches the other transmission band 
edge forming a practically closed trajectory with almost 
zero period of oscillations, see Fig. 0Jc). Note that near 
the left band edge the energy flow in the jg-direction is 
negative but remains positive near the right band edge. 
Further increase of the wavenumber results in the op- 
posite process with the curl forming near the right band 
edge, see Fig.fJJd). For the propagation constants k y cor- 
responding to the second turning point the curl vanishes, 
with a spike forming near the right band edge, and the 
total energy flow along the structure is in the opposite 
direction to the propagation constant. For ky > 2, the 
spike vanishes, and smooth reflection from the right band 
edge is formed. The trajectory again becomes continuous 
and without any ambiguities. 



V. NUMERICAL SIMULATIONS 

For the finite structures with non-adiabatic change of 
the unit-cell size, the approximation of geometric optics 
is not valid, however we can still apply the Bloch-wavc 
formalism for describing the beam propagation, including 
the beam diffraction and formation of interference pat- 
terns. To find the field distribution in the structure, we 
calculate the eigenmodes Ei{z) and eigenvalues k y of a fi- 
nite stack of layers. There are several known approaches 
to solve this problem. The first approach is based on 
the representation of the field in each slab in the form 
of the counter-propagating waves. Boundary conditions 
between the layers define the ratio of the corresponding 
wave amplitudes in the slabs. These amplitudes can be 
represented in terms of either transfer matrix or scat- 
tering matrix formalisms [13, [S|. The transfer matrix 
approach can be used for calculating the eigenmodes and 
eigenvalues for the structures with a small number of pe- 
riods, and it is not suitable for calculating evanescent 
modes in long enough structures. The scattering matrix 
approach provides a better convergence, but it also can- 
not be applied for calculating evanescent modes in the 
structures with more than approximately 100 slabs. An- 
other approach is based on the discretization of the wave 
equation across the structure, and we have used it previ- 
ously to find the eigenvalues of the discrete problem [l5[ . 
However, this approach is computationally intensive, and 
it is not practical for calculating the transmission of very 
long structures. 

Here we implement another approach based on the 
Bloch wave formalism. In Sec. [H]we have already demon- 
strated that for the slowly varying period of the structure 
(i.e., SA <C A) the eigenmodes can be described by Eq. §5§ 
with the amplitudes A\ and wavenumbers k y . We as- 
sume that the structure with a linearly changing period 
is placed between two semi-infinite periodic structures, 
with the period coinciding with that of the adjacent lay- 
ers of the structure, see Fig. [TJ In each unit cell the field 
can be presented as a superposition of the forward and 
backward propagating Bloch waves, 

E%(y,z) = A A U + (z)exp(-i(K b z + k y y)) + (22) 
B A U-(z) exp(-i(-K b z + k y y)), 

where A\ and B\ are the slowly varying amplitudes, U± 
and Kt, are defined by Eqs. $Q and Q, respectively. 

In the left semi-infinite structure, the fields decay ex- 
ponentially, and we choose the amplitudes in the first 
layer as A Al — 1 and B Al — 0. In the right semi- infinite 
structure the forward wave vanishes, and this means that 
A a,. = and B Ar = 1- We calculate the amplitudes in 
the adjacent unit cells using the boundary conditions. 
Repeating this procedure we find the amplitudes in the 
whole structure (for more details, see Appendix A). We 
note that K b is real in the band, hence the Bloch waves 
are propagating, and the modes do not grow exponen- 
tially. In a small number of layers adjacent to the band, 
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FIG. 5: (Color online) Top: ratio between the amplitudes 
of the counter-propagating Bloch waves across the structure 
calculated for k y — 1.838. Bottom: energy flow averaged by 
the unit cell across the structure, for k y = 1.838. 



Kb is imaginary and, consequently, the amplitudes grow 
exponentially. We choose relatively small number of unit 
cells with the parameters corresponding to the bandgap, 
so that the amplitudes and B\ remain reasonable 
across the structure. We note that the absolute value of 
the ratio of the amplitudes of the counter-propagating 
waves in the band is unity, so that the energy is localized 
in the transverse direction inside the structure. We trace 
the dynamics of the Bloch wave amplitudes from left and 
right boundaries separately and compare the amplitudes 
in the middle of the band region, see Fig. \5\ For the 
cigcnmodcs with the eigenvalues k y , the amplitudes of 
the waves at the left and right edges of the structures 
become linearly dependent, i.e., 



det 



1 

Bi/Ai 



1 

By I -Ay 



0. 



(23) 



A=(A r +A,)/2 



Using the procedure described above, we find the propa- 
gation constants k y corresponding to the eigenvalues and 
eigenmodes of the problem. Then, the initial field dis- 
tribution is represented as a superposition of the eigen- 
modes, and the full structure of the fields can be re- 
trieved. 

We calculate the eigenmodes for the structure consist- 
ing of 2100 unit-cells with the width changing linearly 
from the value A; = 4.23 to the value A r = 4.44, and 
SA ~ 2xl0~ 4 . We note that < A >= (A r +A/)/2 ~ 2/3A, 
where A is a free space wavelength. The width of the 
whole structure is approximately 1500A, which is, e.g., 
about few millimeters in infrared regime. We have cho- 
sen the initial field distribution in the form of Gaussian 
beam, see Eq. launched near the right band bound- 
ary of the structure, see Fig. [Sj The width of the beam 
is a = 20(A; + A r )/2, i.c with about 20 unit-cells exited, 
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FIG. 6: (Color online) Distribution of the field amplitude 
in the structure, for the beam with the width a — 20A, at 
k a y = 1.838. Oscillation period is about 23,5 x 100. 



ky = 1.838. We decompose the initial field distribution 
into the superposition of eigenmodes of the structure by 
the least-squares method. The spectrum of the excited 
cigcnmodcs is shown in Fig. GUright). The spectrum is 
narrow, and it is centered near k y = 1.838. The spec- 
trum is to some accuracy equidistant, meaning that the 
field restores its shape after the distance $ = 2tt /5k y . We 
plot the field amplitude averaged over the unit cell. The 
beam curling predicted in the framework of the geometric 
optics is clearly observed. At the beam self-crossing point 
we observe an interference pattern created by the forward 
and backward propagating waves. The field distribution 
is restored almost completely after the first period of os- 
cillations. 

Now we calculate the distribution and direction of the 
Poynting vector. In Fig. [7] we plot the Poynting vector 
averaged across the unit cells. It is clearly seen that near 
the left boundary the direction of the Poynting vector 
along y-axis is opposite to the direction near the right 
boundary. At the self-crossing point corresponding to 
singularity of the group velocity, we observe that the en- 
ergy flow vanishes in the y-dircction. 

We also calculate the average y-component of the 
Pointing vector, defined as < S y >, on the unit cell and 
show its dynamics across the structure for the mode with 
the propagation constant k y = 1.838, see Fig. [SJ It is 
clearly seen that at the left boundary of the band the 
energy flow direction is opposite to that near the right 
boundary, and it vanishes in the center of the structure. 
This dynamics is consistent with the behavior predicted 
within geometric optics approximation in Secs.lITland llVI 
We note that the total energy flow, i.e. the integrated 
flow across the structure, remains positive for this mode 
and for all modes contributing to the beam spectrum. 
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FIG. 7: (Color online) Distribution of the averaged field am- 
plitude and Poynting vector in the chirped structure. Arrows 
show the direction of the energy flow. 



VI. CONCLUSIONS 

We have presented a systematic analysis of the propa- 
gation of electromagnetic waves in chirped periodic struc- 
tures composed of two kinds of alternating layers, the 
layers of negative-index metamaterial and the layers of 
conventional dielectrics, under the condition of the zero 
averaged refractive index. We have considered the chirp 
in the structure parameters introduced by varying the 
thickness of all layers linearly across the structure, and 
we have applied the methods of the geometric optics for 
analyzing the beam propagation and Bloch oscillations 
in such infinite composite structures. For the adiabat- 
ically changing period of structure, we have predicted 
the beam self-crossing in the bands, and have found that 
the energy flow in such multi-layer structures with meta- 
materials may have the opposite direction at the band 
edges, in a sharp contrast to the similar structures com- 
posed of conventional dielectrics. This novel effect of the 
beam curling can be useful for the beam steering in the 
transmission band. 



APPENDIX A 

In this Appendix we discuss the relation between the 
Bloch amplitudes A rn and B m for the m-th layer and the 
amplitudes Aq and Bq in the 0-th layer. 

We describe the field in the m-th unit cell of the pe- 
riodic structure with the period A m as a Bloch wave 
E m {z,y) composed of the counter-propagating compo- 
nents with the amplitudes A m and B m , see Eq. (|22|) . We 
consider the boundary between m-th and (m+l)-th unit 
cells, i.e. z = z m+ i, see Fig. Q] The boundary condition 
at this interface requires that 



Em(z m +i) — E m+ i(z m+ i) 



(Al) 



_1_ dE m (z) 



02 



1 



8E m+1 (z) 



02 



Substituting Eq. |3]) into Eq. (|A1[) . we obtain: 



A m (aj 



B m ( ai + 6, ) = 



= A m+1 (a+e ik - A + b+e~ ik '- A )e- iK >> A + 
+B m+1 (a~e lk ^ A + b~e~ lk " A )e lK " A , 

%[An(,a+ -bt) + B m {ai -fcp] = 



(A2) 



[A„+i(a+e^ A 



b+t 



-ik r K 



)e 



-iK b A 



+ B m+l {a~e ik - A - b-e- lk " A )e zKbA ] 



where ± corresponds to the forward and backward prop- 
agating Bloch waves, respectively. Equation (|A2[) can be 
presented in the matrix form, 



A n 
B n 



= D r 



A m +i 
B m +i 



(A3) 



Using this matrix equation, we can present the field am- 
plitudes at the m-th layer thought the corresponding am- 
plitudes at the 0-th layer, 
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